{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Content under Creative Commons Attribution license CC-BY 4.0, code under BSD 3-Clause License © 2018  by D. Koehn, notebook style sheet by L.A. Barba, N.C. Clementi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Execute this cell to load the notebook's style sheet, then ignore it\n",
    "from IPython.core.display import HTML\n",
    "css_file = '../../style/custom.css'\n",
    "HTML(open(css_file, \"r\").read())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerical stability, dispersion and anisotropy of the 2D acoustic finite difference modelling code\n",
    "\n",
    "Similar to the 1D acoustic FD modelling code, we have to investigate the stability and dispersion of the 2D numerical scheme. Additionally, the numerical dispersion shows in the 2D case an anisotropic behaviour.\n",
    "\n",
    "Let's begin with the CFL-stability criterion ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## CFL-stability criterion for the 2D acoustic FD modelling code\n",
    "\n",
    "As for the 1D code, the maximum size of the timestep $dt$ is limited by the Courant-Friedrichs-Lewy (CFL) criterion:\n",
    "\n",
    "\\begin{equation}\n",
    "dt \\le \\frac{dx}{\\zeta v_{max}}, \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "where $dx$ denotes the spatial grid point distance and $v_{max}$ the maximum P-wave velocity. The factor $\\zeta$ depends on the used FD-operators, dimension and numerical scheme.\n",
    "\n",
    "As for the [1D case](http://nbviewer.jupyter.org/github/daniel-koehn/Theory-of-seismic-waves-II/blob/master/04_FD_stability_dispersion/1_fd_stability_dispersion.ipynb), we estimate the factor $\\zeta$ by the von Neumann analysis, starting with the finite difference approximation of the 2D acoustic wave equation\n",
    "\n",
    "\\begin{equation}\n",
    " \\frac{p_{j,l}^{n+1} - 2 p_{j,l}^n + p_{j,l}^{n-1}}{\\mathrm{d}t^2} \\ = \\ vp_{j,l}^2\\biggl(\\frac{p_{j,l+1}^{n} - 2 p_{j,l}^n + p_{j,l-1}^{n}}{\\mathrm{d}x^2} + \\frac{p_{j+1,l}^{n} - 2 p_{j,l}^n + p_{j-1,l}^{n}}{\\mathrm{d}z^2}\\biggr),\n",
    "\\end{equation}\n",
    "\n",
    "and assuming harmonic plane wave solutions for the pressure wavefield of the form:\n",
    "\n",
    "\\begin{equation}\n",
    "p = exp(i(k_x x + k_z z -\\omega t)),\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "with $i^2=-1$, the wavenumbers $(k_x, k_z)$ in x-/z-direction, respectively, and the circular frequency $\\omega$. Using a regular grid with \n",
    "\n",
    "$dx = dz = dh,$\n",
    "\n",
    "discrete spatial coordinates  \n",
    "\n",
    "$x_j = j dh,$\n",
    "\n",
    "$z_l = l dh,$\n",
    "\n",
    "and times \n",
    "\n",
    "$t_n = n dt.$\n",
    "\n",
    "we can calculate discrete plane wave solutions at the discrete locations and times in eq. (1):  \n",
    "\n",
    "\\begin{align}\n",
    "p_{j,l}^{n+1} &= exp(-i\\omega dt)\\; p_{j,l}^{n},\\\\\n",
    "p_{j,l}^{n-1} &= exp(i\\omega dt)\\; p_{j,l}^{n},\\\\\n",
    "p_{j+1,l}^{n} &= exp(ik_x dh)\\; p_{j,l}^{n},\\\\\n",
    "p_{j-1,l}^{n} &= exp(-ik_x dh)\\; p_{j,l}^{n},\\\\\n",
    "p_{j,l+1}^{n} &= exp(ik_z dh)\\; p_{j,l}^{n},\\\\\n",
    "p_{j,l-1}^{n} &= exp(-ik_z dh)\\; p_{j,l}^{n},\\\\\n",
    "\\end{align}\n",
    "\n",
    "Inserting eqs. (2) - (7) into eq. (1), division by $p_{j,l}^{n}$ and using the definition: \n",
    "\n",
    "\\begin{equation}\n",
    "\\cos(x) = \\frac{exp(ix) + exp(-ix)}{2},\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "yields:\n",
    "\n",
    "\\begin{equation}\n",
    "cos(\\omega dt) - 1 = vp_{j,l}^2 \\frac{dt^2}{dh^2}\\biggl(\\{cos(k_x dh) - 1\\} + \\{cos(k_z dh) - 1\\}\\biggr).\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "Some further rearrangements and division of both sides by 2, leads to:\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{1 - cos(\\omega dt)}{2} = vp_{j,l}^2 \\frac{dt^2}{dh^2}\\biggl(\\frac{1 - cos(k_x dh)}{2} + \\frac{1 - cos(k_z dh)}{2}\\biggr).\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "With the relation \n",
    "\n",
    "\\begin{equation}\n",
    "sin^2\\biggl(\\frac{x}{2}\\biggr) = \\frac{1-cos(x)}{2}, \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "we get \n",
    "\n",
    "\\begin{equation}\n",
    "sin^2\\biggl(\\frac{\\omega dt}{2}\\biggr) = vp_{j,l}^2 \\frac{dt^2}{dh^2}\\biggl(sin^2\\biggl(\\frac{k_x dh}{2}\\biggr)+sin^2\\biggl(\\frac{k_z dh}{2}\\biggr)\\biggr). \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "Taking the square root of both sides finally yields \n",
    "\n",
    "\\begin{equation}\n",
    "sin\\biggl(\\frac{\\omega dt}{2}\\biggr) = vp_{j,l} \\frac{dt}{dh}\\sqrt{sin^2\\biggl(\\frac{k_x dh}{2}\\biggr)+sin^2\\biggl(\\frac{k_z dh}{2}\\biggr)}.\n",
    "\\end{equation}\n",
    "\n",
    "This result implies, that if the Courant number\n",
    "\n",
    "\\begin{equation}\n",
    "\\epsilon = vp_{j,l} \\frac{dt}{dh} \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "is larger than $1/\\sqrt{2}$, you get only imaginary solutions, while the real part is zero (think for a second why). Consequently, the numerical scheme becomes unstable, when the following CFL-criterion is violated\n",
    "\n",
    "\\begin{equation}\n",
    "\\epsilon = vp_{j,l} \\frac{dt}{dh} \\le \\frac{1}{\\sqrt{2}} \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "Rearrangement to the time step dt, assuming that we have defined a spatial grid point distance dh and replacing $vp_{j,l}$ by the maximum P-wave velocity in the FD model $v_{max}$, leads to \n",
    "\n",
    "\\begin{equation}\n",
    "dt  \\le \\frac{dh}{\\sqrt{2}v_{max}}. \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "Therefore, the factor $\\zeta$ in the general CFL-criterion\n",
    "\n",
    "\\begin{equation}\n",
    "dt  \\le \\frac{dh}{\\zeta vp_j}, \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "for the FD solution of the 2D acoustic wave equation using the temporal/spatial 3-point operator to approximate the 2nd derivative is $\\zeta = \\sqrt{2}$. \n",
    "\n",
    "Let's check if this result ist correct:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Import Libraries \n",
    "# ----------------\n",
    "import numpy as np\n",
    "from numba import jit\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from pylab import rcParams\n",
    "\n",
    "# Ignore Warning Messages\n",
    "# -----------------------\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [],
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Definition of modelling parameters\n",
    "# ----------------------------------\n",
    "xmax = 5000.0 # maximum spatial extension of the 1D model in x-direction (m)\n",
    "zmax = xmax   # maximum spatial extension of the 1D model in z-direction (m)\n",
    "dx   = 10.0   # grid point distance in x-direction (m)\n",
    "dz   = dx     # grid point distance in z-direction (m)\n",
    "\n",
    "tmax = 0.8    # maximum recording time of the seismogram (s)\n",
    "dt   = 0.0010 # time step\n",
    "\n",
    "vp0  = 3000.  # P-wave speed in medium (m/s)\n",
    "\n",
    "# acquisition geometry\n",
    "xr = 2000.0   # x-receiver position (m)\n",
    "zr = xr       # z-receiver position (m)\n",
    "\n",
    "xsrc = 2500.0 # x-source position (m)\n",
    "zsrc = xsrc   # z-source position (m)\n",
    "\n",
    "f0   = 20. # dominant frequency of the source (Hz)\n",
    "t0   = 4. / f0 # source time shift (s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# FD_2D_acoustic code with JIT optimization\n",
    "# -----------------------------------------\n",
    "@jit(nopython=True) # use Just-In-Time (JIT) Compilation for C-performance\n",
    "def FD_2D_acoustic_JIT(dt,dx,dz,f0):        \n",
    "    \n",
    "    # define model discretization\n",
    "    # ---------------------------\n",
    "\n",
    "    nx = (int)(xmax/dx) # number of grid points in x-direction\n",
    "    print('nx = ',nx)\n",
    "\n",
    "    nz = (int)(zmax/dz) # number of grid points in x-direction\n",
    "    print('nz = ',nz)\n",
    "\n",
    "    nt = (int)(tmax/dt) # maximum number of time steps            \n",
    "    print('nt = ',nt)\n",
    "\n",
    "    ir = (int)(xr/dx)      # receiver location in grid in x-direction    \n",
    "    jr = (int)(zr/dz)      # receiver location in grid in z-direction\n",
    "\n",
    "    isrc = (int)(xsrc/dx)  # source location in grid in x-direction\n",
    "    jsrc = (int)(zsrc/dz)  # source location in grid in x-direction\n",
    "\n",
    "    # Source time function (Gaussian)\n",
    "    # -------------------------------\n",
    "    src  = np.zeros(nt + 1)\n",
    "    time = np.linspace(0 * dt, nt * dt, nt)\n",
    "\n",
    "    # 1st derivative of a Gaussian\n",
    "    src  = -2. * (time - t0) * (f0 ** 2) * (np.exp(- (f0 ** 2) * (time - t0) ** 2))\n",
    "\n",
    "    # Analytical solution\n",
    "    # -------------------\n",
    "    G    = time * 0.\n",
    "\n",
    "    # Initialize coordinates\n",
    "    # ----------------------\n",
    "    x    = np.arange(nx)\n",
    "    x    = x * dx       # coordinates in x-direction (m)\n",
    "\n",
    "    z    = np.arange(nz)\n",
    "    z    = z * dz       # coordinates in z-direction (m)\n",
    "\n",
    "    # calculate source-receiver distance\n",
    "    r = np.sqrt((x[ir] - x[isrc])**2 + (z[jr] - z[jsrc])**2)\n",
    "\n",
    "    for it in range(nt): # Calculate Green's function (Heaviside function)\n",
    "        if (time[it] - r / vp0) >= 0:\n",
    "            G[it] = 1. / (2 * np.pi * vp0**2) * (1. / np.sqrt(time[it]**2 - (r/vp0)**2))\n",
    "    Gc   = np.convolve(G, src * dt)\n",
    "    Gc   = Gc[0:nt]\n",
    "\n",
    "    # Initialize model (assume homogeneous model)\n",
    "    # -------------------------------------------\n",
    "    vp    = np.zeros((nx,nz))\n",
    "    vp2    = np.zeros((nx,nz))\n",
    "\n",
    "    vp  = vp + vp0       # initialize wave velocity in model\n",
    "    vp2 = vp**2\n",
    "    \n",
    "    # Initialize empty pressure arrays\n",
    "    # --------------------------------\n",
    "    p    = np.zeros((nx,nz)) # p at time n (now)\n",
    "    pold = np.zeros((nx,nz)) # p at time n-1 (past)\n",
    "    pnew = np.zeros((nx,nz)) # p at time n+1 (present)\n",
    "    d2px = np.zeros((nx,nz)) # 2nd spatial x-derivative of p\n",
    "    d2pz = np.zeros((nx,nz)) # 2nd spatial z-derivative of p\n",
    "\n",
    "    # Initialize empty seismogram\n",
    "    # ---------------------------\n",
    "    seis = np.zeros(nt) \n",
    "    \n",
    "    # Calculate Partial Derivatives\n",
    "    # -----------------------------\n",
    "    for it in range(nt):\n",
    "    \n",
    "        # FD approximation of spatial derivative by 3 point operator\n",
    "        for i in range(1, nx - 1):\n",
    "            for j in range(1, nz - 1):\n",
    "                \n",
    "                d2px[i,j] = (p[i + 1,j] - 2 * p[i,j] + p[i - 1,j]) / dx**2                \n",
    "                d2pz[i,j] = (p[i,j + 1] - 2 * p[i,j] + p[i,j - 1]) / dz**2\n",
    "\n",
    "        # Time Extrapolation\n",
    "        # ------------------\n",
    "        pnew = 2 * p - pold + vp2 * dt**2 * (d2px + d2pz)\n",
    "\n",
    "        # Add Source Term at isrc\n",
    "        # -----------------------\n",
    "        # Absolute pressure w.r.t analytical solution\n",
    "        pnew[isrc,jsrc] = pnew[isrc,jsrc] + src[it] / (dx * dz) * dt ** 2\n",
    "                \n",
    "        # Remap Time Levels\n",
    "        # -----------------\n",
    "        pold, p = p, pnew\n",
    "    \n",
    "        # Output of Seismogram\n",
    "        # -----------------\n",
    "        seis[it] = p[ir,jr]        \n",
    "        \n",
    "    return time, seis, Gc, p     # return last pressure wave field snapshot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To seperate modelling and visualization of the results, we introduce the following plotting function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compare FD Seismogram with analytical solution\n",
    "# ---------------------------------------------- \n",
    "def plot_seis(time,seis_FD,seis_analy):     \n",
    "    \n",
    "    # Define figure size\n",
    "    rcParams['figure.figsize'] = 12, 5\n",
    "    \n",
    "    plt.plot(time, seis_FD, 'b-',lw=3,label=\"FD solution\") # plot FD seismogram\n",
    "    plt.plot(time, seis_analy,'r--',lw=3,label=\"Analytical solution\") # plot analytical solution\n",
    "    \n",
    "    plt.xlim(time[0], time[-1])\n",
    "    \n",
    "    plt.title('Seismogram')\n",
    "    plt.xlabel('Time (s)')\n",
    "    plt.ylabel('Amplitude')\n",
    "    plt.legend()\n",
    "    plt.grid()\n",
    "    \n",
    "    plt.show() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dx   = 10.0   # grid point distance in x-direction (m)\n",
    "dz   = dx     # grid point distance in z-direction (m)\n",
    "f0   = 20     # centre frequency of the source wavelet (Hz)\n",
    "\n",
    "# define \n",
    "zeta = np.sqrt(2)\n",
    "\n",
    "# calculate dt according to the CFL-criterion\n",
    "dt = dx / (zeta * vp0)\n",
    "\n",
    "%time time, seis_FD, seis_analy, p = FD_2D_acoustic_JIT(dt,dx,dz,f0)\n",
    "plot_seis(time,seis_FD,seis_analy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical Grid Dispersion\n",
    "\n",
    "While the FD solution above is stable, it is subject to some numerical dispersion when compared with the analytical solution. The grid point distance $dx = 10\\; m$, P-wave velocity $vp = 3000\\; m/s$ and a maximum frequency $f_{max} \\approx 2*f_0 = 40\\; Hz$ leads to ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fmax = 2 * f0 \n",
    "N_lam = vp0 / (dx * fmax)\n",
    "print(\"N_lam = \",N_lam)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$N_\\lambda = 7.5$ gridpoints per minimum wavelength. Let's increase it to $N_\\lambda = 12$, which yields ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_lam = 12 \n",
    "dx   = vp0 / (N_lam * fmax)\n",
    "dz   = dx     # grid point distance in z-direction (m)\n",
    "f0   = 20     # centre frequency of the source wavelet (Hz)\n",
    "\n",
    "# define \n",
    "zeta = np.sqrt(2)\n",
    "\n",
    "# calculate dt according to the CFL-criterion\n",
    "dt = dx / (zeta * vp0)\n",
    "\n",
    "%time time, seis_FD, seis_analy, p = FD_2D_acoustic_JIT(dt,dx,dz,f0)\n",
    "plot_seis(time,seis_FD,seis_analy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... an improved fit of the 2D analytical by the FD solution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical Anisotropy\n",
    "\n",
    "Compared to the 1D acoustic case, the numerical dispersion behaves a little bit differently in the 2D FD approximation. To illustrate this problem, we model the pressure wavefield for $t_{max} = 0.8\\; s$for a fixed grid point distance of $dx = 10\\;m$ and a centre frequency of the source wavelet $f0 = 100\\; Hz$ which corresonds to $N_\\lambda = 1.5$ grid points per minimum wavelength."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define dx/dz and calculate dt according to the CFL-criterion\n",
    "dx   = 10.0   # grid point distance in x-direction (m)\n",
    "dz   = dx     # grid point distance in z-direction (m)\n",
    "# define zeta for the CFL criterion\n",
    "zeta = np.sqrt(2)\n",
    "dt = dx / (zeta * vp0)\n",
    "\n",
    "f0 = 100\n",
    "time, seis_FD, seis_analy, p = FD_2D_acoustic_JIT(dt,dx,dz,f0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot last pressure wavefield snapshot at Tmax = 0.8 s\n",
    "# -----------------------------------------------------\n",
    "rcParams['figure.figsize'] = 8, 8 # define figure size\n",
    "\n",
    "clip = 1e-7  # image clipping\n",
    "extent = [0.0, xmax/1000, 0.0, zmax/1000]\n",
    "\n",
    "# Plot wavefield snapshot at tmax = 0.8 s\n",
    "plt.imshow(p.T,interpolation='none',cmap='seismic',vmin=-clip,vmax=clip,extent=extent)\n",
    "plt.title('Numerical anisotropy')\n",
    "plt.xlabel('x (km)')\n",
    "plt.ylabel('z (km)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clearly, a spatial discretization with $N_\\lambda = 1.5$ gridpoints per minimum wavelength leads to modelling results significantly affected by dispersion. However, note that the dispersion is direction dependent. Maximum dispersion coincides with the main Cartesian coordinate axes, while the dispersion reaches a minimum at 45° angles relative to the main axes. This problem is called **Numerical Anisotropy** and can be explained by the angle dependence of the numerical phase velocity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What we learned:\n",
    "\n",
    "* Estimation of the CFL-criterion for the FD approximation of the 2D acoustic wave equation using von Neumann analysis\n",
    "* Verification of CFL and grid dispersion criterion by 2D numerical FD modelling\n",
    "* Be aware of the influence of numerical anisotropy on FD modelling results"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
